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1  Introduction 

This  report  documents  progress  in  the  implementation  of  two  new  types  of  source  functions  in 
the  WAVEWATCH  III®  model.  The  WAVEWATCH  model  was  originally  developed  at  Delft 
University  (Tolman  1991).  In  its  current  form,  it  is  referred  to  as  WAVEWATCH  III®  (“WW3”), 
developed  at  NOAA’s  NCEP  (Tolman  et  al.  2002).  At  time  of  writing,  the  last  public  release  of 
WW3  was  WW3  version  3.14  (Tolman  2009),  and  WW3  version  4  is  under  active  development 
via  a  Subversion  (svn)  software  versioning  and  revision  control  system  administered  by  NCEP. 

The  governing  equation  of  WW3  and  most  other  “third  generation”  or  “3G”  wave  models  is  the 
action  balance  equation.  Simplified  here  from  the  WW3  form  for  purposes  of  presentation,  the 
action  balance  equation  is: 


dt 


dC  N  dC  N  dCN  dCsN 
+  — - —  +  — —  + — —  +  -  6 


dx 


dy 


da 


de 


S_ 

<7 


(1) 


The  prognostic  variable  is  the  wave  action  density  N,  equal  to  energy  density  divided  by  angular 
relative  frequency  (N  —  E /o'),  and  is  a  function  of  space  and  time,  N  =  N{x,y,9,a,t).  Relative 
frequency  a  is  the  wave  frequency  measured  from  a  frame  of  reference  moving  with  a  current,  if 
a  current  exists;  9  is  wave  direction;  C  is  the  wave  action  propagation  speed  in  (x,  y,  9,  o,  t) 
space.  In  absence  of  currents,  Cx  is  the  x-component  of  the  group  velocity  C„.  The  right  hand  side 
of  the  governing  equation  is  the  total  of  source/sink  terms  expressed  as  rate  of  change  of  wave 
action  density,  where  S  =  S(x,y,  9,  a,  t )  is  most  generally  represented  by  three  terms,  =  Sltl  + 
$nl  +  Sds:  input  by  wind,  nonlinear  interactions,  and  dissipation,  respectively. 

Sea  ice  and  mud  affect  the  length  and  dissipation  rate  of  wind- generated  ocean  waves.  The  ice- 
modified  (or  mud-modified)  wavenumber  can  be  expressed  as  a  complex  number  k  =  kr  +  ikt, 
with  the  real  part  kr  representing  impact  of  the  sea  ice/mud  on  the  physical  wavelength  and 
propagation  speeds,  producing  something  analogous  to  shoaling  and  refraction  by  bathymetry, 
whereas  the  imaginary  part  of  the  complex  wavenumber,  kt,  is  an  exponential  decay  coefficient 
ki(x,  y,  t,  a)  (depending  on  location,  time  and  frequency,  respectively),  representing  wave 
attenuation,  and  can  be  introduced  in  a  wave  model  such  as  WW3  as  Sice/E  —  —2Cgki  ,  where 
Sice  is  one  of  several  dissipation  mechanisms,  along  with  whitecapping,  for  example,  Sds  — 

Swc  +  Sice  +  Smud  +  •■•  on  the  right-hand  side  of  the  governing  equation  (see  also  Komen  et  al. 
(1994,  pg.  170)).  The  kr  modified  by  ice/mud  would  enter  the  model  via  the  C  calculations  on 
the  left-hand  side  of  the  governing  equation.  Though  the  procedure  is  non-trivial  with  regard  to 
necessary  code  changes,  especially  with  regard  to  i/o,  the  fundamentals  are  straightforward,  e.g. 
Rogers  and  Holland  (2009  and  subsequent  unpublished  work)  modified  a  similar  model,  SWAN 
(Booij  et  al.  1999)  to  include  the  effects  of  a  viscous  mud  layer  using  the  same  approach 
(k  —  kr  +  ikt)  previously.  The  wa ve-mud  interaction  theory  implemented  by  Rogers  and 
Holland  (2009)  was  derived  by  Ng  (2000)  for  a  viscous  mud  layer  under  a  (very  weakly)  viscous 
water  layer.  This  “two  layer”  approach  is  taken  in  some  solutions  for  wave-ice  interaction  theory, 
e.g.  Keller  (1998),  though  he  assumes  that  the  water  layer  is  inviscid,  and  of  course  the  second 
layer  is  above  rather  than  below  the  water  layer. 
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In  this  report,  we  describe  the  implementation  of  new  routines  to  represent  the  effect  of  ice  and 
mud  on  waves  in  the  WAVEWATCH  III  model.  Both  implementations  utilize  the  concept  of 
complex  wavenumber  k  —  kr  +  ikt.  However,  only  the  imaginary  component  is  addressed  in 
WW3  in  this  report  (thus  it  is  more  limited  than  what  was  done  for  mud  in  the  SWAN  model). 
Modification  of  the  real  part  has,  at  time  of  writing,  not  been  addressed  yet  in  WW3,  though  this 
is  part  of  our  plans. 

General  code  modifications  are  described  in  Section  2.  Rather  than  provide  further  background 
information  for  the  two  source  terms  here,  we  introduce  them  in  their  respective  sections.  Section 
3  deals  with  the  effect  of  ice  on  waves.  Section  4  deals  with  the  effect  of  mud  on  waves. 
Suggestions  for  further  work  are  discussed  in  Section  5. 

2  WW3  code:  input  methods 

With  regard  to  model  coding,  the  most  challenging  task  associated  with  this  project  so  far  has 
been  not  in  the  source  term  routines  themselves,  but  rather  in  the  code  associated  with  the 
processing  of  user  input.  The  latter  is  necessary  since  the  new  source  terms  require  new  variables 
to  be  input  by  the  user.  In  the  case  of  mud,  we  introduce  new  variables:  mud  thickness,  mud 
viscosity,  and  mud  density.  In  the  case  of  ice,  we  allow  up  to  five  new  parameters.  These  can  be 
referred  to  generically  as  Cice  l,  Cice2,  ...,  Cice  5.  In  the  code,  Cice  l  is  referred  to  as  “ICECOEF1” 
(a  local  scalar  parameter  in  the  source  function  routine)  or  “ICEP1”  (a  2D  array  shared  to  the 
source  function  routine  via  a  module  which  describes  the  spatial  variation  of  the  parameter  on 
the  computational  grid).  Our  intent  here  is  to  allow  the  meaning  of  the  ice  parameters  to  vary 
depending  on  which  Sice  routine  is  selected.  For  example,  in  an  implemented  Sice  routine 
(described  below,  referred  to  as  the  “Liu  routine”),  Cice  l  represents  the  ice  thickness  (in  meters) 
and  Cice2  represents  the  eddy  viscosity  in  the  turbulent  boundary  layer  beneath  the  ice. 

Some  remarks  about  this  strategy: 

1)  External  variables  already  available,  like  currents,  water  depth,  wind,  ice  concentration 
can  also  be  used  (though  probably  only  ice  concentration  is  useful  for  Sice)  and  do  not 
count  against  the  maximum  total  of  5  parameters. 

2)  External  variables  not  already  available,  like  water  temperature,  salinity,  ice  thickness, 
effective  viscosity,  could  be  used  but  would  count  against  the  maximum  total  of  5 
parameters. 

3)  If  a  developer  feels  that  the  five  ice  parameters  and  three  mud  parameters  are  insufficient, 
these  numbers  can  be  increased,  but  unfortunately,  this  is  not  easy  to  code.  We  point  out 
that  since  these  are  rather  specific  and  specialized  physics  routines,  there  probably  will 
not  be  many  situations  in  which  ice  and  mud  are  used  in  the  same  simulation.  A 
developer  can  potentially  exploit  this  by  using  mud  “parameter  space”  for  ice,  as  Cice  6, 
^ice, C[ce  s- 

4)  Ideas  for  potential  use  of  the  ice  parameters  are  discussed  in  Section  5. 

The  new  parameters  are  read  in  using  the  same  methods  that  already  exist  for  other  scalar 
parameters,  such  as  ice  concentration  and  water  level.  The  variables  are  allowed  to  vary  in  time 
and  space.  In  case  of  spatially  varying  parameters,  these  are  read  in  via  the  ww3_prep  program, 
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using  instructions  in  the  ww3_prep.inp  user-input  file.  Otherwise,  the  user  can  use  the  simpler 
option  of  specifying  them  as  homogeneous  (but  potentially  time-varying)  fields  via  the 
ww3_shel  program,  using  instructions  in  the  ww3_shel.inp  user-input  file.  Though  the  latter 
method  is  simpler,  it  is  not  expected  to  find  much  use  other  than  for  idealized  test  cases.  The 
ww3_prep  approach  of  WW3  supports  a  number  of  different  methods  of  user  input.  For 
example,  the  user  can  provide  the  ice  parameters  as  ascii  files  on  a  non-WW3  grid,  and 
ww3_prep  will  interpolate  in  time  and  space  to  the  WW3  computational  grid(s). 

The  WW3  code  and  test  cases  described  in  this  report  are  kept  on  the  NRL  svn 
repository,  which  was  last  synchronized  with  the  trunk  of  the  NCEP  svn  repository  at  revision 
21198  (Sep.  19  2012).  The  latest  revision  to  the  NRL  svn  repository  at  time  of  writing  was  228. 

3  Effect  of  ice  on  waves 

The  implementation  of  Sice  is  described  in  this  section. 

3.1  Background:  wave  modeling  in  the  Arctic 

The  mutual  interactions  between  ocean  waves  and  sea  ice  cover  play  a  crucial  role  for  planning 
safe  operations  in  the  Arctic  Ocean.  Therefore,  wave  and  ice  interactions  should  be  among  the 
center-pieces  of  the  operational  wave  forecasting  system.  A  research  objective  of  NRL  and  the 
Office  of  Naval  Research  (ONR)  is  to  study  these  interactions  in  the  marginal  ice  zone  (MIZ)  of 
the  Arctic  Ocean,  and  develop  techniques  for  modeling  the  effect  of  sea  ice  on  wave  energy  and 
wavelength.  A  number  of  theories  and  models  have  been  developed  to  describe  this 
phenomenon,  e.g.  Keller  (1998),  Liu  et  al.  (1991),  Squire  et  al.  (1995),  Wadhams  et  al.  (1986), 
Wadhams  et  al.  (1988),  and  Wang  and  Shen  (2010).  A  brief  review  of  these  methods  is  given  in 
the  Appendix  (Section  8). 

The  retreating  ice  cover  implies  an  increase  in  fetch  for  generation  of  waves  in  the  Arctic.  This, 
combined  with  more  frequent  incoming  cyclones  in  the  Arctic  (Sepp  and  Jaagus  2011)  naturally 
leads  to  more  severe  wave  conditions.  The  reduction  of  the  permanent  polar  pack  ice  also  implies 
that  regions  of  the  Arctic  that  could  previously  be  ignored  in  operational  numerical  wave  models 
must  now  be  considered.  Lortunately,  NRL  has  recently  extended  the  capability  of  the 
WAVEWATCH  III  (WW3,  Tolman  1991,  2009)  model  so  that  it  can  be  applied  on  irregular 
grids  (Rogers  and  Campbell  2009).  An  implementation  of  WW3  for  the  Arctic  has  been 
successfully  demonstrated  in  the  beta  queue  at  LNMOC  (bleet  Numerical  Meteorology  and 
Oceanography  Center),  on  the  same  grid  as  used  for  the  Arctic  atmospheric  model  (COAMPS, 
Hodur  1997),  with  significantly  better  resolution  (15-20  km) — and  better  forcing — than  that  of 
the  global  wave  model.  Ligure  1  and  Ligure  2  below  shows  examples  of  test  simulations  on  this 
grid.  During  LY 11,  NRL  has  extended  WW3  to  allow  use  of  two-way  nesting  “mosaic”  approach 
of  Tolman  (2008)  with  curvilinear  grids.  Realtime  surface  current  and  ice  concentration  values 
are  available  from  the  1/12°  Arctic  Cap  Nowcast/Lorecast  System  (ACNLS),  developed  at  NRL 
(Posey  et  al.  2010).  Lurther,  within  the  next  two  years,  funded  by  the  Earth  System  Prediction 
Capability  (ESPC)  program,  NRL  will  begin  transitioning  to  the  Naval  Oceanographic  Office 
(NAVOCEANO)  an  Earth  System  Model  Lramework  (ESML)-based,  coupled  wave-ocean 
model  on  a  global  high  resolution  tripolar  grid.  Both  wave  model  grid  methods,  the  curvilinear 
two-way  nested  Arctic  regional  grid  and  the  global  tripolar  (curvilinear)  grid,  address  the 
traditional  problems  associated  with  extending  a  regular  global  grid  to  high  latitudes,  e.g.  the 


3 


operational  WW3  at  FNMOC  stops  at  78°  latitude.  With  all  this  in  mind,  we  can  see  that  from  a 
technical  standpoint,  the  operational  Navy  is  well-positioned  for  forecasting  waves  in  the  Arctic 
over  the  next  decade. 


Figure  1.  Propagation  test  with  WAVEWATCH  III  model  on  curvilinear  Arctic  grid.  [Significant 
waveheight  Hs  in  meters]  No  ice,  winds,  or  boundary  forcing  are  included,  and  the  region  above 
89°  is  treated  as  land  to  avoid  directional  singularity.  The  initial  condition  (geographic 
distribution)  is  a  Gaussian  spike  in  the  wave  field.  The  plotted  condition  is  after  several  hours  of 
propagation. 
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Figure  2.  Source  term  test  with  WAVEWATCH  III  model  with  global  and  Arctic  grid.  The 
curvilinear  Arctic  grid  is  shown  here.  [Significant  waveheight  Hs  in  meters]  Ice,  winds,  and 
boundary  forcing  are  included  in  this  hindcast.  This  is  the  result  at  May  25  2009  12Z,  after  a  12 
hour  simulation  (from  cold  start).  In  this  simulation,  two-way  nesting  is  performed,  such  that 
wave  spectra  from  this  nest  can  propagate  across  the  boundary  into  the  global  model,  and  vice 
versa. 

Unfortunately,  the  situation  with  regard  to  the  physics  of  wave  models  in  the  Arctic  is  much  less 
optimistic.  The  key  physical  process,  wave  attenuation  by  interaction  with  ice  floes  in  the 
Marginal  Ice  Zone  (MIZ),  is  treated  within  the  propagation  routine  of  the  model,  with  the  percent 
transmission  of  wave  energy  through  ice  being  a  simple  function  of  ice  concentration.  There  is 
no  connection  to  any  physical  mechanism  for  wave  attenuation;  this  artificial  “dissipation”  is  not 
dependent  on  frequency  and  has  a  spurious  dependence  on  grid  aspect  ratio  (Section  3.2).  This 
simple,  non-physical  approach  could  nevertheless  be  justified  on  the  grounds  that  operational 
characterization  of  the  ice  is  limited  (ice  concentration  only)  and  further  that  the  existing 
physical  mechanisms  available  from  the  literature  which  might  be  implemented  in  the  forecast 
model  are  1)  too  numerous  and  too  varied  to  select  from  and  2)  too  poorly  informed,  e.g.  what  is 
the  value  of  a  theoretical  model  which  requires  input  parameters  that  are  impossible  to  estimate? 
Compounding  the  problem  is  a  very  limited  number  of  studies  estimating  attenuation  from 
observations,  which  would  normally  be  used  to  calibrate  and  verify  a  numerical  model. 

As  mentioned  above,  a  specific  NRL  research  objective  is  to  develop  techniques  for  modeling 
the  effect  of  sea  ice  on  waves.  In  the  present  effort,  we  utilize  modeling  codes  that  are  currently 
used  operationally — WW3  in  the  case  of  the  wave  model — distinguishing  our  aim  from  more 
detailed  process-based  modeling  investigations,  such  as  models  of  individual  waves  and  ice 
floes. 


5 


3.2  Background:  existing  methods  vs.  new  methods  of  representing  the 
effect  of  ice  on  waves 

The  existing  method  of  treating  ice  is  described  by  Tolman  (2003,  2009).  Ice  concentration 
e(x,  y,  t)  is  specified  by  the  user,  in  addition  to  two  constant  parameters  which  describes  the 
minimum  concentration  which  affects  the  waves,  ec  0  and  the  concentration  at  which  wave 
energy  is  completely  blocked,  ecn.  For  example,  the  parameters  might  be  set  at  ec0  —  0.25  and 
ec,n= 0-75.  For  ice  concentrations  between  ec  0  and  ecn,  the  wave  energy  is  partially 
blocked/transmitted  based  on  linear  interpolation  between  the  two  values:  the  cell  transparency 
in  the  x  direction  is  calculated  as  ax(x,  y,  t)  —  (ec  n  —  e(x,  y,  t)) / {ec,n  ~  ec,o)-  This  method 
necessitates  grid-specific  calibrations;  this  characteristic  is  perhaps  best  illustrated  by 
recognizing  that  as  Ax  0  ,  the  method  will  give  infinite  dissipation. 

In  the  v3.14  public  release  of  WW3  (Tolman  2009),  the  amount  of  blocking  had  an 
unfortunate  additional  dependence  on  grid  aspect  ratio,  but  this  was  removed  in  the  development 
version  (v4),  by  Dr.  Ardhuin  (Ifremer).  A  contemporary  change  was  to  add  an  option  to  replace 
this  ax  calculation  with  a  new  one:  ax(x,y,  t )  =  exp  (— e(x,  y,  t)Ax)/Le,  where  Le  (denoted 
“FICEL”  and  “LICE”  in  the  code)  is  a  new  user-specified  constant  parameter,  apparently  a 
dissipation  length  scale.  It  can  be  shown  that  for  e=l  and  Le  =  1/kj  this  method  is  a  numerical 
approximation  of  the  analytical  formula  given  by  (5)  below.  This  further  implies  that  the  method 
is  effectively  similar  to  our  first  method  described  in  Section  3.3.1,  but  using  variable  ice 
concentration  e(x,y,  t)  and  the  additional  constant  parameter  Le  in  place  of  the  variable 
parameter  kj(x,y,  t).  The  original  parameters  ec  0  and  ec  n,  are  not  used  by  this  new,  optional 
method.  This  method  is  not  documented. 

In  any  case,  these  existing  methods  represent  the  dissipation  of  waves  by  interaction  with 
sea  ice  using  the  LHS  of  (1)  (propagation/blocking),  and  our  objective  here  is  to  move  this  to  the 
RHS  of  (1)  (dynamics). 

Reflection  by  icebergs  (distinct  from  sea  ice)  is  added  to  the  development  version  of  the 
code  by  Dr.  F.  Ardhuin.  This  new  feature  is  documented  in  the  development  version  of  the 
manual.  This  is  added  primarily  to  address  the  overprediction  of  significant  waveheight  (SWH) 
by  global  operational  models  near  Antarctica. 

3.3  WW3  code:  method  of  including  Sice 

The  methods  of  user-input  has  already  been  explained  above.  In  this  section,  we  describe  two 
methods  of  calculating  Sice  which  we  have  implemented  in  WW3.  It  is  anticipated  that  additional 
methods,  e.g.  Keller  (1998),  Wang  and  Shen  (2010),  will  be  implemented  in  the  future. 

In  comparison  to  the  Tolman  (2003)  method  of  representing  the  dissipation  by  ice  as  a 
per-cell  partial  blocking  mechanism  on  the  LHS  of  equation  (1),  the  new  method  is  to  treat  as  a 
source  term  on  the  RHS  of  equation  (1).  The  new  method  has  the  advantage  of  removing  the 
proportional  dependence  of  dissipation  on  resolution.  As  Ax  -»  0  ,  the  old  method  would  give 
infinite  dissipation,  whereas  the  new  method  converges  to  a  proper,  continuous  solution. 

Technical  details.  To  follow  the  WW3  convention,  each  Sice  method  would  have  a 
Fortran  file  associated  with  it.  However,  to  simplify  code  during  the  development  process,  the 
two  Sice  methods  are  kept  in  the  same  Fortran  file  (w3siclmd.ftn)  for  now,  and  the  user  selects 
the  Sice  method  using  a  namelist  variable.  At  a  later  stage,  these  will  be  expanded  into  multiple 
Fortran  files  (w3siclmd.ftn,  w3sic2md.ftn  w3sic3md.ftn,  etc.)  which  are  selected  via  “switch” 
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file,  following  WW3  convention.  Though  the  latter  approach  unfortunately  tends  to  result  in 
substantial  repeated/redundant  code,  it  is  highly  beneficial  when  multiple  groups  develop  source 
term  methods. 

3.3.1  Method  1:  dissipation  rate  constant  in  frequency  space. 

The  first  implemented  method  is  for  the  user  to  specify  y,  t)  which  is  uniform  in  frequency 
space,  C[ce  l  =  /Cj.  In  this  case,  the  amount  of  information  read  in  has  not  changed  from  the 
Tolman  (2003)  method  of  using  ice  concentration,  e(x,y,  t). 

3.3.2  Method  2:  Liu  et  al.  (1991) 

This  method  is  based  on  the  papers  by  Liu  and  Mollo-Christensen  (1988)  and  Liu  et  al.  (1991); 
these  will  be  denoted  as  “LMC”  and  “LHV”  here.  This  is  a  model  for  “viscous  attenuation  by  a 
sea  ice  cover”,  derived  on  the  assumption  that  dissipation  is  caused  by  turbulence  in  the 
boundary  layer  between  the  ice  floes  and  the  water  layer,  with  the  ice  modeled  as  a  continuous 
thin  elastic  plate.  As  mentioned  above,  input  ice  parameters  are  ice  thickness  (in  meters)  and  an 
“eddy  viscosity  in  the  turbulent  boundary  layer  beneath  the  ice”,  v.  Cice  l  represents  the  former 
and  Cice  2  represents  the  latter.  Ice  concentration  is  not  an  input  to  this  routine;  this  is  discussed 
further  below. 

A  description  of  the  code  follows: 

1.  General  routine.  If  non-zero  Cice  l,  the  forward  dispersion  routine  (item  2  below)  is 
called.  This  is  used  to  calculate  a,  the  spatial  exponential  decay  rate  of  energy1.  From 
there,  D  —  —2Cgk^  and  finally  Sice  —  DN.  Here,  D  represents  the  temporal  decay  rate, 

D  =  Sice/E .  Recall  from  above  that  Cg  is  group  velocity,  S  is  the  source  term  (following 
WAMDI  (1988)  convention),  E  is  spectral  energy  density  and  N  =  E /a  is  spectral  action 
density.  The  variable  D  varies  in  frequency  space  but  is  constant  in  directional  space.  The 
variable  Sice  varies  in  directional  space  via  dependence  on  N.  Further,  recall  that  Cice 
parameters  vary  in  geographic  space  and  time,  e.g.  Cice  l{x,y,  t). 

2.  Forward  dispersion  relation  routine.  This  is  very  much  like  a  traditional  dispersion 
relation:  given  frequency  /,  find  wavenumber  k.  However,  in  this  case,  the  wavenumber 
is  a  complex  number.  The  LHV  dispersion  relation  cannot  be  solved  directly,  so  Newton- 
Raphson  method  is  used  here,  calling  the  “reverse  dispersion”  routine  (item  3  below), 
which  is  directly  solvable.  Inputs  to  this  routine  are  (all  being  local  terms):  ice  thickness, 
eddy  viscosity,  water  depth,  and  /.  Output  from  this  routine  are:  kr ,  Cg,  a.  The  first  two 
can  potentially  be  used  later  to  feed  back  to  model  kinematics,  to  produce  refraction  (in 
case  of  kr)  and  shoaling  (in  case  of  Cg)  by  ice.  Only  a  is  relevant  to  the  Sice  calculation. 

3.  Reverse  dispersion  routine.  This  is  a  directly  solvable  calculation  using  the  equations  of 
LHV.  Inputs  to  this  routine  are:  ice  thickness,  eddy  viscosity,  water  depth,  kr.  Outputs 
from  this  routine  are:  /,  Cg ,  and  a. 

Equations: 

The  reader  is  referred  to  LHV,  equations  on  page  4606.  The  key  equations  are: 


1  Note:  a  is  exponential  decay  rate  of  energy  while  /q  is  exponential  decay  rate  of  amplitude,  so  lcL  =  a/2.  There  is 
no  intended  connection  between  a  and  ax ,  the  latter  being  the  variable  used  by  Tolman  (2009)  to  represent  cell 
transparency. 
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a2  =  ( gkr  +  B/c^)/(coth(/cr/iw)  +  krM ) 


(2) 


Cg  =  (g  +  (  5  +  4/crM)B/c^)/(2cr(l  +  krM )2) 


(3) 


a:  =  (Vvor/cr)/(<25V 2(1  +  krM )) 


(4) 


In  our  notation,  hw  is  water  depth  and  h,  is  ice  thickness.  There  is  an  apparent  typo  in  equation 
(1)  of  LHV,  coth  (kh{)  should  be  coth  ( khw ).  The  variables  B  and  M  quantify  the  effects  of  the 
bending  of  the  ice  and  inertia  of  the  ice,  respectively.  Both  of  these  variables  depend  on  hL  (for 
these  equations,  see  LMC). 

Example  calculations  of  dissipation  rate  a  using  the  LHV  model  are  shown  in  Figure  3.  In  this 
case,  the  three  described  routines  are  coded  in  Matlab,  but  they  have  been  re-coded  in  Fortran  for 
the  purpose  of  application  in  WW3. 


x  10  6 


Figure  3.  Example  calculations  of  dissipation  rate  a  using  the  LHV  model,  coded  using  the  three 
routines  as  described  in  the  text.  Units  of  hw  and  are  in  meters;  units  of  v  are  m"/sec. 
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LHV  state:  “The  only  tuning  parameter  is  the  turbulent  eddy  viscosity,  and  it  is  a  function  of  the 
flow  conditions  in  the  turbulent  boundary  layer  which  are  determined  by  the  ice  thickness,  floe 
sizes,  ice  concentration,  and  wavelength.”  The  eddy  viscosity  term  v  given  by  LHV  is 
unfortunately  “highly  variable”  (their  words),  and  “not  a  physical  parameter”,  which  suggests 
that  it  is  difficult  to  specify  in  practice.  In  LHV,  many  values  are  referenced  and  used2: 

1.  v=160.0x  10-4  m2/sec  (Brennecke  1921) 

2.  v=24.0x  1(T4  m2/sec  (Hunkins  1966) 

3.  v=3450.0x  1CT4  m2/sec  (LHV  Fig.  11) 

4.  v=4.0x  1(T4  m2/sec  (LHV  Fig.  12) 

5.  v=150.0x  10“4  m2/sec  (LHV  Fig.  13) 

6.  v=54.0x  1(T4  m2/sec  (LHV  Fig.  14) 

7.  v=384.0x  10“4  m2/sec  (LHV  Fig.  15) 

8.  v=1536.0x  1(T4  m2/sec  (LHV  Fig.  16) 

Another  criticism  of  this  source  term  is  that  it  does  not  use  the  ice  concentration  in  actual 
calculations.  The  model  assumes  a  continuous  ice  layer  (100%  concentration),  so  the  method 
appears  to  simply  rely  on  concentration  being  high:  “When  the  ice  is  highly  compact  with  high 
concentration,  the  flexural  waves  obey  the  dispersion  relation. .  .as  similar  waves  in  a  continuous 
ice  sheet.”  Later,  “Five  of  these  cases  with  high  ice  concentration  (larger  than  60%)  in  the  MIZ 
have  been  selected”.  For  general  use,  it  would  be  better  to  include  concentration  in  the 
calculations.  This  might  be  added  by  incorporating  concentration  as  a  scaling  factor. 

Other  settings;  all  three  are  from  LHV,  pg.  4606  are: 

1)  Young's  modulus  of  elasticity  is  set  to  E  =  6.0  x  109  N/m2 

2)  Poisson's  ratio  is  set  to  s  =  0.3. 

3)  The  relation  between  ice  density  and  water  density:  pL  =  0.9pw. 

3.4  One-dimensional  tests  with  Sice 

Figure  4  shows  results  for  a  simple  one-dimensional  test  case.  The  analytical  expression  used 
here  is: 

tf(x)  =  H0e-ktx 

where  H  is  significant  waveheight  and  H0  is  significant  waveheight  at  x  —  0.  In  Figure  4, 

Ax  =  1  km  is  used.  The  model  captures  the  decay  well  for  weaker  dissipation  values,  but  at  the 
highest  dissipation  values,  the  model  decay  is  somewhat  slower  than  the  analytical  solution. 

Even  with  the  highest  dissipation  rate,  the  error  might  be  considered  tolerable.  To  demonstrate 
sensitivity  to  geographic  resolution,  results  with  Ax  =  10  km  are  shown  in  Figure  5.  In  this  case, 
the  numerical  error  for  the  higher  dissipation  rates  is  clearly  not  acceptable.  These  results  suggest 
that  expected  dissipation  rates  must  be  part  of  the  decision  with  respect  to  what  resolution  to  use: 
if  is  large,  then  the  spatial  resolution  cannot  be  too  coarse,  or  the  numerical  representation  is 
poor. 


2  Note:  In  our  implementation,  the  user  specifies  eddy  viscosity  in  units  of  m2/sec  even  though  values  are  given  in 
units  of  cm2/sec  in  LHV. 
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dx  =  1000m 


Figure  4.  One-dimensional  tests  using  Ax  =  1  km.  Significant  waveheight  is  plotted.  Waves  are 
initialized  at  the  left  boundary,  x  =  0.  Results  using  various  settings  are  shown.  Units  of 
are  1/m.  The  dissipation  rate  kL  is  stationary  and  constant  in  x  and  in  frequency  space.  Solid 
lines:  calculations  using  an  analytical  expression.  Circles:  WW3  output.  This  test  case  is 
included  in  the  NRL  svn  repository  for  WW3. 
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dx  =  10km 


Figure  5.  Identical  to  Figure  4,  except  that  Ax  =  10  km  is  used. 

The  question  of  “how  much  error  is  tolerable”  is,  of  course,  subjective.  If  we  assume  that  2.7% 
error  in  Hm0  is  intolerable  in  our  idealized  simulations,  experiments  with  a  number  of  resolutions 
can  be  summarized  as  follows  (with  kt  given  in  m"1). 

•  for  Ax=20.0  km,  kt  should  not  exceed  3.5e-6 

•  for  Ax=5.0  km,  kt  should  not  exceed  2.0e-5 

•  for  Ax=2.5  km,  kt  should  not  exceed  5.0e-5 

•  for  Ax=1.0  km,  /q  should  not  exceed  2.0e-4  (assumes  2.7%  Hs  error  is  intolerable) 

•  for  Ax=0.35  km,  error  is  less  than  2.1%  for  all  kt  tested  (up  to  1.0e-3) 

•  for  Ax=0.10  km,  error  is  less  than  1.3%  for  all  kt  tested  (up  to  1.0e-3) 

For  reference,  Arctic  Cap  Nowcast/Forecast  System  (Posey  et  al.  2010)  is  1/12°,  so  Ax=9.25  km 
north-south. 

The  analytical  formula  also  allows  us  to  put  the  dissipation  rate  of  the  original  model  (partial 
blocking  due  to  ice  concentration)  in  context  with  kt  values.  For  example,  if  the  original  model 
has  a  grid  cell  with  50%  blocking,  and  has  resolution  of  55  km,  this  is  equivalent  to  kt  —  1.26  x 
10"5  m1. 
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3.5  Two-dimensional  tests  with  Sice 

Two  two-dimensional  test  cases  have  been  added  to  the  svn  repository  for  WW3.  The  first,  like 
the  one-dimensional  test  case  of  the  previous  section,  uses  a  /q  value  read  in  from  input  files 
which  is  constant  in  frequency  space.  The  domain  is  square,  with  nx  —  ny  =  51  and  Ax  —  Ay  — 
5  km.  Boundary  forcing  is  uniform  and  steady  along  the  southern  and  eastern  boundaries, 
producing  swells  propagating  from  the  southeast  ( 9  =  135°).  The  initial  condition  is  a  uniform 
wave  field  equivalent  to  the  boundary  forcing.  In  the  northwest  quadrant  of  the  domain,  ice 
appears  and  disappears  during  the  duration,  with  /q  =  0,  /q  —  1  x  10-5  m'1,  /q  =  2  x  10~5  m  , 
kt  —  1  x  10-5  m and  then  /q  =  0.  Results  from  this  test  case  are  not  shown,  but  are  consistent 
with  expectations:  the  non- stationary,  non-uniform  ice  specification  is  validated  to  work  as 
intended. 

The  second  two-dimensional  test  case  uses  the  Liu  et  al.  (1991)  attenuation  model,  with 
ice  thickness  of  1  m  and  ice  eddy  viscosity  parameter  ve  =  1  x  10-4  m2/sec.  The  computational 
grid  uses  three  frequencies:  0.08  Hz,  0.10  Hz,  and  1.25  Hz,  with  most  of  the  energy  in  the 
specified  boundary  forcing  corresponding  to  the  central  frequency.  With  this  model,  /q  is,  of 
course,  frequency-dependent.  For  these  three  frequencies  and  ice  parameters,  /q  =  6.0  x  10~6, 
8.9  x  10~6,  and  9.5x  10-6  m  1  respectively.  The  computational  grid  is  square,  with  nx  —  101, 
ny  =  51,  Ax  =  5  km,  and  Ay  =  10  km.  As  in  the  other  two-dimensional  test,  the  boundary 
forcing  is  uniform  and  steady  along  the  southern  and  eastern  boundaries,  producing  swells 
propagating  from  the  southeast  (9  —  135°).  The  initial  condition  consists  of  near-zero  wave 
energy  on  the  interior  of  the  domain.  Ice  is  again  specified  in  the  northwest  quadrant,  but  unlike 
the  other  two-dimensional  test,  in  this  test,  ice  conditions  are  stationary.  Swell  propagates  into 
the  ice  and  eventually  a  steady  state  is  reached,  with  total  simulation  duration  being  30  hours. 
The  final  state  is  shown  in  Figure  6. 
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Figure  6.  Final  state  of  the  second  two-dimensional  test  case,  as  explained  in  the  text. 

4  Effect  of  mud  on  waves 

Wave  damping  by  muddy  seabeds  is  generally  understood  to  occur  when  wave-generated 
stresses  exceed  the  limiting  strength  of  the  bed,  causing  the  liquefaction  of  some  or  all  of  the 
mud  layer  into  a  viscous  “fluid  mud”.  Internal  waves  are  then  generated  at  the  water-mud 
interface,  and  their  energy  is  dissipated  relatively  rapidly  by  viscosity  within  the  fluid  mud  layer. 
Several  theoretical  approaches  have  been  used  to  model  the  effect  of  mud  on  waves,  representing 
the  mud  as  a  purely  viscous  fluid  (e.g.,  Dalrymple  and  Liu,  1978;  Ng,  2000;  Winterwerp  et  al., 
2007),  or  alternatively  as  viscoelastic  (e.g.,  MacPherson,  1980;  Jiang  and  Mehta,  1996;  Zhang 
and  Ng,  2006)  or  plastic  (e.g.,  Mei  and  Liu,  1987). 

Until  recently,  the  effect  of  mud  on  waves  has  had  little  or  no  representation  in  generally 
available  nearshore  wave  models  such  as  SWAN  (Booij  et  al.,  1999)  and  WW3,  which  have  been 
limited  to  rigid  beds.  This  has  been  a  major  shortcoming,  forcing  users  to  specify  unrealistic 
bottom  friction  parameters  in  muddy  areas  in  order  to  get  the  desired  wave  dissipation 
characteristics. 

The  present  section  describes  the  implementation  in  WW3  of  two  viscous  fluid  mud  dissipation 
formulations,  borrowing  heavily  from  very  similar  implementations  added  to  the  SWAN  wave 
model  by  Rogers  and  Holland  (2009). 


13 


4.1  Background:  implementation  in  SWAN 

Three  formulations  for  wave  damping  by  purely  viscous  mud  have  been  implemented  in  SWAN. 
The  first  formulation  (Rogers  and  Holland,  2009,  based  on  Dalrymple  and  Liu,  1978)  treats  the 
mud  as  a  laminar  viscous  fluid.  This  is  a  relatively  accurate  method,  but  it  requires  a  complex 
iterative  technique  which  significantly  lengthens  computational  time.  The  second  formulation 
(Rogers  and  Holland,  2009,  based  on  Ng,  2000)  is  simplified  relative  to  Dalrymple  and  Liu 
(1978)  in  that  it  assumes  the  mud  layer  to  be  thin.  This  formulation  computes  a  mud- induced 
group  velocity  from  the  real  part  of  the  mud-induced  wavenumber,  which  allows  the  effects  of 
mud  on  wave  refraction,  shoaling,  and  de-shoaling  to  be  estimated.  The  third  formulation 
(Winterwerp  et  al.,  2007)  integrates  the  energy  transport  across  the  water/mud  interface  over  one 
wave  period,  based  on  earlier  work  by  Gade  (1958)  and  De  Wit  (1995).  Unlike  the  other 
implementations,  it  assumes  the  water  to  be  inviscid  and  does  not  consider  mud  effects  on  wave 
phase  and  group  velocities. 

4.2  WW3  code:  method  of  including  Smud 

For  the  present  project,  two  of  the  above  formulations  for  the  dissipation  of  wave  energy  by 
viscous  mud  were  implemented  in  WW3.  The  first  implementation  (module  “w3sbt8md”)  is 
based  on  the  formulation  of  Dalrymple  and  Liu  (1978).  A  second  implementation  (module 
“w3sbt9md”)  is  based  on  the  formulation  of  Ng  (2000).  The  code  that  was  originally  created  for 
these  formulations  in  the  spectral  wave  model  SWAN  was  transferred,  with  a  small  number  of 
modifications,  directly  into  the  WW3  modules.  Additional  modifications  were  made  to  several 
other  WW3  subroutines  to  allow  users  to  turn  on/off  the  mud  dissipation  routines  and  to  input 
field  data  for  mud  thickness,  density,  and  kinematic  viscosity  (Section  2). 

A  description  of  the  WW3  code  follows: 

1.  General  routines.  If  the  user  specifies  “BT8”  or  “BT9”  in  the  “switch”  parameter  file,  the 
preprocessor  will  activate  code  statements  calling  the  w3sbt8  or  w3sbt9  subroutines  in 
the  computation  of  source  terms  by  module  w3srcemd.  For  non-uniform  input  fields,  the 
user  must  create  “ww3_prep”  input  files  for  each  of  the  mud-related  parameters:  mudd 
(density),  mudt  (thickness),  and  mudv  (viscosity).  Field  data  on  the  distribution  of  mud 
parameter  values  throughout  the  grid  are  read  either  from  these  files  or  from  separate 
field  parameter  files  referenced  in  the  prep  files. 

2.  Dalrymple  &  Liu  routine  (w3sbt8).  The  wave  dissipation  by  the  fluid  mud  is  computed 
using  an  iterative  procedure  that  converges  to  the  complex  mud-induced  wave  number, 
kmud .  Dissipation  due  to  mud  at  each  frequency  is  determined  from  the  imaginary  part  of 
this  wave  number  as 

D„u,d  =  2  ■  imag(kmud )  •  Cg  mui  ^ 

where  C  mud  is  the  mud-induced  wave  group  velocity.  This  dissipation  is  added  to 

contributions  from  other  source/sink  terms  by  w3srce.  For  additional  details,  see 
Dalrymple  &  Liu  (1978). 
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3. 


Ng  routine  (w3sbt9).  After  initialization  of  mud  field  data  and  assignment/computation  of 
various  local  parameters,  the  w3sbt9  routine  calls  the  “Ng”  subroutine  to  compute  the 
mud-induced  dissipation.  The  Ng  routine  determines  two  second-order  coefficients 
( Br ,  Bt )  that  are  then  used  to  compute  wave  attenuation  due  to  mud  as  follows: 


Dmud  =  imag(k2)  = 


sinh  IkJ'i  +  2k  Jt 


(7) 


where  Sm  is  the  Stokes  boundary  layer  thickness  for  mud,  h  is  water  depth,  and  k,  and 


k2  are  the  first-  and  second-order  parts  of  the  wave  number,  respectively,  in  a  Taylor 
expansion  about  the  mud- water  interface.  For  additional  details,  see  Ng  (2000).  Results 
from  the  Ng  subroutine  are  returned  to  w3sbt9,  which  passes  them  back  to  w3srce  where 
the  mud- induced  dissipation  is  again  added  to  contributions  from  other  source/sink  terms. 


Note  a  fixed  value  for  water  kinematic  viscosity  “nu_water”  (1.31E-6  nr/s)  was  added  to  the 
module  “constants”  (i.e.,  file  “constants. ftn”).  This  is  the  kinematic  viscosity  of  pure  water  at 
10°C  and  is  in  accordance  with  the  water  density  value  (1000  kg/m3)  that  is  used  throughout 
WW3.  This  value  for  kinematic  viscosity  of  water  is  used  in  both  w3sbt8  and  w3sbt9. 


4.3  One-dimensional  test  cases  with  Smud 

Test  cases  (“mud_testl”  and  “mud_test2”,  respectively)  were  created  for  WW3  to  simulate  ID 
wave  propagation  for  a  distance  of  100  km  over  a  flat  bottom  with  a  mud  layer  of  constant 
thickness  using  each  of  the  above  formulations.  Parameters  for  both  tests  are  generally  set  to 
match  those  used  in  Fig.  3  of  Rogers  and  Holland  (2009),  except  for  mud  thickness  and  grid 
size/spacing.  These  quantities  are  allowed  to  vary  in  order  to  investigate  their  effects  on  model 
accuracy. 

The  following  is  a  list  of  features  of  both  ID  test  cases: 

•  spectral,  spatial,  and  time  settings: 

o  three  (3)  frequencies  from  0.08  to  0.125  Hz;  24  directions 
o  initial  wave  height  =  1 .0  m 
o  nx=24-120,  ny=3 
o  Ax =Ay=  1.0  to  5.0  km 
o  boundary  forcing  from  west  boundary 
o  boundary  forcing:  9=270°  (waves  from  southeast), fp=0. 10  Hz 
o  Stalling  time  :  1968/06/06  00:00:00  UTC 
o  Ending  time  :  1968/06/06  12:00:00  UTC 

•  mud  parameters  are  constant  for  entire  domain: 

o  mud  density=  1310  kg/m 
o  mud  thickness=0.01  -  0.4  m 
o  mud  viscosity=7.60E-03  rrr/s 

The  objective  of  these  tests  is  to  compare  the  actual  decay  calculated  by  WW3  with  the  methods 
of  Dalrymple  and  Liu  (1978)  and  Ng  (2000)  to  the  expected  exponential  decay.  For  the 
comparison,  multiple  values  were  used  for  mud  thickness,  which  caused  the  the  exponential 
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decay  coefficient,  k,.  to  range  from  0.23e-05  to  12.74e-05  m  As  seen  earlier  with  ice 
dissipation,  model  accuracy  was  affected  by  the  grid  spacing  used  in  these  calculations.  Figure  8 
shows  estimated  H,„o  values  (circles  and  asterisks)  plotted  versus  theory  (i.e.,  Eq.  5;  solid  lines) 
for  grid  spacing  of  Ax  =  1  km,  while  Figure  8  shows  the  same  comparisons  for  grid  spacing  of 
Ax  =  5  km. 


j= 

o 

E 

X 


dx  =  1km 
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Figure  7.  Comparison  of  WW3-predicted  significant  wave  heights  (circles  and  asterisks)  with 
expected  exponential  decay  (solid  lines)  for  several  values  of  C,  using  grid  spacing  of  Ax  =  1  km 
and  mud  parameter  settings  described  in  the  text. 


In  general,  the  error  levels  for  the  two  methods  are  relatively  similar.  To  limit  the  maximum 
error  for  these  simulations  to  less  than  2.7%  (to  be  consistent  with  Section  3.4),  the  value  of  k; 
should  not  exceed  roughly  4.5e-05  m  1  (for  1-km  spacing)  or  1.25e-05  m  1  (for  5-km  spacing). 

It  is  noted  that  the  computational  time  required  by  the  Dalrymple  and  Fiu  (1978)  formulation  - 
roughly  30  seconds  -  is  twenty  times  greater  than  that  required  by  the  Ng  (2000)  formulation  - 
roughly  1.5  seconds. 
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Figure  8.  Same  as  Figure  7,  but  for  Ax  =  5km  spacing. 


4.4  Two-dimensional  test  case  with  Smud 

This  case  is  similar  to  the  2d  Sjce  test  case.  The  following  is  a  “bullet  list”  of  features  of  this  test 
case: 

•  features  that  are  the  same  as  S,-ce  test: 

o  three  (3)  frequencies  from  0.08  to  0.125  Hz 
o  nx=101,ny=51 
o  Ax=5  km,  Ay=10  km 

o  boundary  forcing  from  south  and  east  boundaries 
o  boundary  forcing:  0=135°, ^;=0. 10  Hz 
o  Stalling  time  :  1968/06/06  00:00:00  UTC 
o  Ending  time  :  1968/06/07  06:00:00  UTC 

•  patch  of  mud  in  northwest  quadrant  of  domain  (x=0  to  25  km,  y=25  to  50  km) 

•  mud  density=  1310  kg/m 

•  mud  thickness=0.4  m 

•  mud  viscosity=7.60E-03  m2/s 
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The  two-dimensional  test  case  displayed  a  similar  difference  in  required  computational  time  for 
the  two  methods  to  that  seen  earlier  in  ID.  For  these  simulations,  the  w3sbt8md  module  based 
on  Dalrymple  and  Liu  (1978)  required  9  minutes  39  seconds  on  a  workstation  using  three  cpu 
cores.  In  contrast,  the  w3sbt9md  module  based  on  Ng  (2000)  required  only  15.1  seconds.  Thus, 
in  this  case  the  Ng  (2000)  formulation  was  roughly  38  times  faster  than  that  of  Dalrymple  and 
Liu  (1978). 

5  Discussion 

Dissemination  of  code 

WW3  has  always  been  an  open-source  model.  During  the  past  several  years,  the  model  has 
started  a  transition  from  a  code  predominantly  authored  by  a  single  person  to  a  “community 
model”.  In  fact,  in  terms  of  current  development,  it  is  unambiguously  a  community  model.  The 
logistics  of  the  collaboration  are  handled  by  a  Subversion  repository  at  NOAA/NCEP.  At 
present,  the  code  described  in  this  report  is  not  on  that  svn  repository,  but  is  instead  on  an  NRL 
repository.  Our  intent  is  to  create  a  new  branch  on  the  NOAA/NCEP  repository  with  this  code 
(and  associated  test  cases)  immediately  after  publication  of  this  report. 

Alternate  theoretical  models  for  Siro. 
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From  a  physical  basis,  the  LHV  model  is  perhaps  not  the  most  credible  model  that  can  be  found 
in  the  literature.  It  was  selected  as  the  first  model  to  implement  because  of  its  relative  simplicity. 
Other  models,  for  example,  require  the  determination  of  multiple  complex  roots,  from  which  one 
must  be  selected.  A  review  of  other  possible  methods  is  given  in  Section  8.  One  or  more  of  these 
will  be  implemented  in  WW3  in  the  future. 

Alternatives  to  calculating  Siro  within  WW3 

As  we  look  at  more  complex  theoretical  models  for  Sice,  we  may  find  that  the  computational  cost 
is  not  practical.  In  particular,  it  may  be  wasteful  to  make  a  large  number  of  highly  similar  yet 
computationally  expensive  calculations  repeatedly:  “highly  similar”  since  in  many  cases,  the  ice 
conditions  will  be  similar  in  neighboring  grid  cells,  or  in  consecutive  time  steps.  If  the  source 
term  is  linear,  then  there  is  no  dependence  of  D  —  S/E  on  spectral  energy  level,  so  the  spatial 
and  temporal  variation  of  wave  conditions  becomes  irrelevant.  Thus,  it  can  be  reasoned,  that  the 
dissipation  rate  D  might  be  calculated  on  a  reduced  grid,  at  reduced  intervals,  and  then  ingested 
into  WW3  for  quick  translation  onto  the  computational  grid.  One  paradigm  is  the  “look-up  table” 
approach  where  WW3  pulls  D  from  a  table  of  pre-calculated  values  based  on  many  possible  ice 
conditions.  Another  paradigm  is  to  provide  WW3  with  5  to  8  parameters,  as  Cice  l  (x,  y,  t), 

Cice  2(x,  y,  t),  . . .,  Cice  8(x,  y,  t),  also  based  on  some  pre-calculations  outside  of  WW3.  Then,  for 
a  given  (x,  y,  t),  the  parameters  Cice  l,  Cice  2,  . . .,  Cice  8  (up  to  8  numbers)  are  used  to  describe  the 
one-dimensional  function  D(f ). 3  The  existing  method  of  i/o  is  more  conducive  to  the  second 
approach. 

The  simplest  approach  would  be  to  read  in  D  (/)  on  a  discrete,  coarse  frequency  grid  (8 
frequencies)  and  interpolate  these  numbers  onto  the  model  computational  grid’s  rif  frequencies 
(typically  between  25  and  40  frequencies  are  used).  Or,  some  fitting  could  be  used,  e.g.  D(/)  = 
a0  +  ci\f  +  a2/2  +  a3/3  +  a4/4.  In  fact,  the  fitting  can  be  to  any  parametric  form  that  can  be 
imagined.  Matlab  is  particular  powerful  in  this  regal'd,  since  it  will  fit  to  a  user- specified 
parametric  form  via  least  squares. 

Finally,  we  note  that  WW3  apparently  does  have  the  capability  to  ingest  data  assimilation 
data  on  the  computational  frequency  grid.  At  time  of  writing,  we  are  not  aware  of  how  mature 
this  piece  of  code  is,  but  potentially  the  same  methods  could  be  used  to  ingest  D  (/)  on  the 
model’s  computational  frequency  grid. 

Future  improvement:  kr 

As  mentioned  above,  we  plan  to  include  the  effect  of  ice  and  mud  on  the  real  pail  of  the 
wavenumber,  kr,  by  passing  the  phase  velocity  and  group  velocity  variables  back  to  the  main 
routines  for  calculation  of  refraction  and  shoaling  effects. 

Unresolved  complication:  other  source  terms  in  presence  of  ice 

It  stands  to  reason  that  the  deepwater  source  terms  Sin  ,  Sds  and  Snl4  do  not  behave  the  same 
under  partial  ice  cover  as  they  do  in  open  water.  The  most  simple  way  to  address  this  would  be  to 
multiply  each  term  by  the  open  water  fraction,  (1  —  e).  This  will  be  our  initial  approach4. 


3  Keep  in  mind  that  D(/),  fed/)  and  a(f  )  are  all  dissipation  rates  and  are  directly  related,  so  if  one  is  known,  the 
others  are  known. 

4  A  time  of  writing,  this  is  not  done  yet,  though  apparently  the  ST4  source  term  package  already  reduces  Sin  in  this 
manner.  Ideally,  this  operation  should  be  performed  at  a  higher  level  than  the  individual  source  term  packages,  to 
ensure  uniformity. 
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However,  this  simple  approach  is  not  necessarily  the  best:  this  warrants  further  consideration  and 
study.  Observational  data  and  more  process-based  modeling  might  help  with  this  question. 

Applicability  of  Srm.n  within  a  conditionally  stable  model 

WW3,  ran  on  structured  grids  (either  regular  or  irregular  types),  is  conditionally  stable.  As  such, 
there  is  a  practical  limit  on  how  fine  the  geographic  resolution  can  be  set.  Often,  this  will 
discourage  applications  in  the  littoral,  such  as  the  application  of  SWAN  by  Rogers  and  Holland 
(2009),  where  Ax  =  252  m  (two-dimensional  case)  and  Ax  =  50  m  (one-dimensional  case), 
where  the  unconditionally  stable — and  optionally  stationary — SWAN  must  be  considered  as  a 
more  efficient  alternative  to  WW3.  This,  in  turn,  makes  physical  representation  of  wave-seafloor 
interaction  a  lower  priority  in  WW3  than  in  SWAN.  Even  so,  these  interactions  can  be 
significant  on  broad  continental  shelves  modeled  at  resolutions  more  appropriate  to  WW3,  e.g. 
Ax  =  4  km;  a  specific  example  would  be  the  North  Carolina  shelf,  which  is  0(80  km)  across, 
e.g.  see  Ardhuin  et  al.  (2003).  Further,  methods  are  being  introduced  to  enable  higher  resolution 
applications  of  WW3.  Van  der  Westhuysen  and  Tolman  (2012)  introduced  a  quasi- stationary 
method  of  computation  that  can  be  regarded  as  a  compromise  (or  hybrid)  between  the  stationary 
mode  of  computation  optionally  used  with  SWAN  and  traditional  time-stepping,  which  was  the 
only  method  available  in  WW3  prior.  Further,  the  unstructured-grid  methods  of  Roland  et  al. 
(2009)  have  been  implemented  in  an  experimental  version  of  WW3.  Both  developments  suggest 
that,  in  the  near'  future,  it  will  be  possible  to  make  computationally  efficient  applications  of 
WW3  at  high  resolution. 

Closing  remarks 

On  the  subject  of  the  Sice  implementation  in  WW3,  and  the  question  “why  was  this  not  done 
sooner?”,  as  noted  earlier,  the  argument  has  been  that  the  treatment  of  ice  as  a  partial  blocking 
mechanism  via  ice  concentration  (Tolman  2003)  is  sensible,  as:  a)  ice  concentration  is  the  only 
operationally  available  ice  variable,  b)  the  method  is  relatively  cheap,  and  c)  the  end  result  is 
very  similar'  to  what  the  more  sophisticated  Sice  might  provide.  We  accept  (a)  and  (b)  as 
legitimate  arguments,  but  (c)  is  doubtful,  since  the  original  method  is  dependent  on  grid 
resolution  and  aspect  ratio,  requiring  grid-specific  recalibration.  Even  (c)  becomes  a  reasonable 
argument  if  one  considers  the  recent  improvements  to  this  pari  of  the  code  by  Ifremer  (see 
Section  3.2).  In  this  regard,  the  Sice  implementation  herein  might  be  regarded  as  somewhat  in 
advance  of  what  would  be  feasible  operationally.  This  is  a  traditional  problem  with  wave 
modeling:  features  become  available  which  will  not  always  be  applicable  due  to  limits  on  input. 
For  example,  a  high-resolution  two-dimensional  littoral  wave  model  will  not  generally  have 
adequate  bathymetry  in  an  operational  context;  this  is  the  traditional  argument  for  using  simple 
one-dimensional  “surf  models”  operationally.  Similar'  statements  could  be  made  regarding  the 
availability  of  mud  thickness,  viscosity,  and  density.  These  are  valid  arguments,  but  two  counter¬ 
arguments  can  be  made.  First,  we  can  anticipate  that  as  technology  external  to  the  wave  model 
improves,  more  information  about  the  ice  will  become  available.  For  example,  ice  thickness  or 
floe  size  distribution  might  be  derived  from  ice  model  or  remote  sensing  in  the  near'  future;  NRF 
is  now  actively  researching  methods  to  derive  ice  thickness  from  satellite.  Second,  these  new 
routines  can  be  applied  in  non-operational  contexts,  where  the  ice  (or  mud)  is  more  adequately 
prescribed,  such  as  the  case  with  hindcasts  corresponding  to  field  experiments,  already  planned 
with  newly  funded  Office  of  Naval  Research  initiatives.  Fastly,  we  point  out  that  the  simpler 
methods,  even  after  the  Ifremer  updates,  do  not  provide  for  dependence  of  attenuation  on  wave 
frequency  or  the  dependence  of  phase  velocity  and  group  velocity  on  ice.  Though  the 
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quantitative  dependencies  still  need  to  be  worked  out,  the  existence  of  such  dependencies  is  clear 
enough  and  should  be  accommodated  in  subsequent  research  codes. 
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8  Appendix:  Brief  review  of  literature  regarding  Sice  calculation  methods 

There  is  no  shortage  of  theoretical  models  for  representing  the  effects  of  waves  on  ice.  The  larger 
challenges  are  to  select  the  most  appropriate  models,  to  apply  them  only  in  applications 
consistent  with  their  underlying  assumptions,  and  to  determine  the  most  suitable  inputs  for  these 
models,  since  they  unfortunately  are  universally  not  framed  in  terms  of  variables  that  are 
available  operationally.  We  give  a  brief,  partial  survey  of  available  models  here.  Though  we  are 
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primarily  interested  in  the  MIZ,  models  for  the  continuous  ice  problem  are  discussed  here.  A 
more  detailed  review  can  be  found  in  Squire  et  al.  (1995),  since  updated  by  Squire  (2007). 


There  are  a  few  broad  categories  of  models  for  the  MIZ.  The  first  that  we  will  discuss  is  the  type 
of  model  based  on  scattering  from  individual  floes.  Squire  et  al.  (1995)  argue  quite  reasonably 
that  this  type  of  model  is  most  appropriately  used  in  situations  of  less  compact  ice.  This  method 
is  given  by  Wadhams  (1973a)  and  described  in  Wadhams  et  al.  (1986,  1988).  This  model 
requires  a  distribution  of  floe  diameter  di  and  thickness  hL.  Masson  and  LeBlond  (1989)  also 
derive  a  “scattering  by  floes”  model  assuming  cylindrical  floes  (circular  from  top  view) 
organized  in  a  hexagonal  pattern.  It  requires  floe  draft,  radius,  water  depth,  and  distance  to  other 
floes.  Perrie  and  Hu  (1996)  based  their  model  on  this  one.  Kohout  and  Meylan  (2008)  have  a 
similar  model,  requiring  the  mean  size  of  floes,  and  ice  concentration.  Dumont  et  al.  (201 1)  use 
this  model. 

The  second  type  is  more  appropriately  used  in  situations  where  the  ice  floes  are  more  compact, 
potentially  colliding.  Here,  the  floes  are  treated  as  a  single  layer  with  specific  rheology5.  Weber 
(1987)  and  Keller  (1998)  treat  this  ice  cover  as  a  viscous  fluid  layer  in  models  for  broken  ice 
(frazil,  brash,  pancake);  they  require  effective  viscosity  (or  eddy  viscosity)  and  density.  Newyear 
and  Martin  (1999)  use  this  model  in  a  study  of  grease  ice.  The  Wang  and  Shen  (2010)  “unified 
rheological  model”  is  similar  to  Keller  (1998)  but  adds  elasticity;  elasticity  must  be  quantified 
via  the  effective  shear  modulus  G.  A  particular  challenge  in  this  type  of  model  is  to  recreate  the 
so-called  “roll-over”  of  dissipation  seen  in  many  observations,  which  is  technical  jargon  for  the 
non-monotonic  dependence  of  the  dissipation  rate  on  wave  period. 

The  third  type  of  model  is  similar  to  the  second,  insofar  as  it  is  intended  for  highly  compact  ice 
(e.g.  shore-fast  ice),  but  instead  of  presenting  the  problem  in  terms  of  two  “fluid”  layers  with  the 
dissipation  being  caused  by  the  rheology  of  the  ice  layer,  here  the  dissipation  is  attributed  to 
“eddy  viscosity  in  the  turbulent  boundary  layer  beneath  the  ice”.  An  example  of  this  model  is  Liu 
and  Mollo-Christensen  (1988)  which  is  used  by  Liu  et  al.  (1991)  and  Liu  et  al.  (1993)  in 
combination  with  a  model  for  continuous  ice,  Wadhams  (1973b). 

The  fourth  type  of  model  for  the  MIZ  is  described  by  Squire  et  al.  (1995)  as  a  “mass-loading 
model”.  The  Wadhams  and  Holt  (1991)  model  is  one  example.  This  type  of  model  is  of  less 
interest  to  the  present  study,  since  it  does  not  predict  attenuation,  i.e.  dissipation  rate  kt  —  0  (see 
Section  1).  However,  the  propagation  velocity  is  affected  via  modification  of  the  real  paid  of  the 
wavenumber  kr,  so  it  can  block  slower  waves  in  a  manner  similar  to  blocking  of  short  wind 
waves  by  opposing  currents. 

The  continuous  ice  models  are  most  often  associated  with  wave  propagation  through  the  central 
Arctic,  though  there  is  applicability  to  shore-fast  ice  and  perhaps  highly  compact  sea  ice. 
Wadhams  (1973b)  presents  one  such  model,  based  on  plastic-elastic  deformation  of  the  ice  sheet 
by  swell;  plastic  deformation  or  “creep”  provides  the  actual  dissipation.  This  is  applied  by  Liu  et 
al.  (1991)  and  Liu  et  al.  (1993)  to  the  MIZ  problem  along  with  the  boundary  layer  model  as 
mentioned  above.  Here  solution  for  k  requires  E  (Young’s  modulus  of  elasticity),  B  (coefficient 


5  We  acknowledge  the  apparent  contradiction  of  treating  discontinuous  ice  as  a  continuous  layer.  Since  the  situation 
being  modeled  is  not  continuous,  we  prefer  not  to  group  this  under  the  category  of  “continuous  ice  model’’. 
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for  bending)  and  M  (coefficient  for  inertia  of  the  ice).  Squire  (1984)  also  presents  a  model  for 
“shore-fast  ice”;  here,  there  is  a  linear  viscoelastic  plate  over  a  fluid;  the  viscous  quality  is 
intended  to  represent  creep  in  the  ice  layer. 

A  second  broad  category  of  continuous  ice  model  is  that  used  by  Vaughan  et  al.  (2009)  and 
Squire  et  al.  (2009).  This  model  is  based  on  scattering  from  discontinuities  in  ice  thickness,  thus 
requiring  ice  thickness  profile  along  axis  of  wave  propagation.  Squire  et  al.  (2009)  use  this 
model  to  predict  dissipation  rates  using  ice  thickness  profiles  from  upward-looking  sonar 
collected  by  a  submarine. 
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